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Abstract 

This paper is devoted to show the results obtained by using the magnitude- 
squared coherence for determining order-chaos transition in a system de- 
scribed by the logistic equation dynamics. For determining the power spectral 
density of a chaotic finite-duration discrete-time sequence the Welch average 
periodogram method was used. This method has the advantage that can 
be applied to any stationary signal by using the discrete Fourier transform 
(DFT) representation of a discrete-time series which allows an effective com- 
putation via fast Fourier transform (FFT) algorithm, and that can be applied 
to a discrete-time series shorter than that required by nonlinear dynamical 
analysis methods. The estimate of the Inverse Average Magnitude-Squared 
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Coherence Index (lAMSCI) for each discrete-time series in the set obtained 
from the logistic mapping was calculated. The control parameter, r, ranges 
in the interval [2.8,4] for producing the discrete-time series set. When the 
condition r > 3.57 is satisfied, each discrete-time series exhibited a positive 
value for lAMSCI estimate, indicating a high level of coherence loss of the 
signal and corresponding to a chaotic behavior of the dynamical system. Its 
effectiveness was demonstrated by comparing the results with those obtained 
by calculating the largest Lyapounov exponent of the time series set obtained 
from the logistic equation. 
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I. INTRODUCTION 



Nonlinear dynamical methods for analyzing a discrete-time series have been widely used 
in several fields of scientific research. For characterizing a scalar discrete-time series metric 
tools of nonlinear dynamical methods, such as dimensions, exponents, and entropies are 
calculated. Applying these tools to embed a discrete-time series in a phase space to obtain a 
phase-space portrait of the attr actor of the dynamical system is required. In order to make 
the discrete-time series embedding defining the two embedding parameters which are the 
embedding dimension, D^, and the time delay, r, is required. Several methods have been 
proposed for determining adequate values of these parameters [4,22,13,10,7,15,11]. When the 
dimensionality, m, of the dynamical system is known, it is directly considered that De ~ m, 
but it is not the case when we only know a scalar discrete-time series obtained by measuring 
one of the accessible variables of the state space of a dynamical system. The embedding 
dimension is usually estimate in accordance with Takens theorem, i.e., > 2 * D2, where 
D2 is the correlation dimension of the attractor in the reconstructed phase space and it is 
required to be also estimate before Dg be known. After determining the two embedding pa- 
rameters, De and r, one proceeds to calculate estimates of the metric tools. Two important 
requirements must be considered at this point. One of these is related to the stationarity of 
the discrete-time series and the other one is related to the length of the data point sequence. 
Nonlinear dynamical methods applying requires that a discrete-time series to be analyzed 
be a stationary free-noise sequence during the observation time. For calculating the correla- 
tion dimension, D2, and the largest Lyapounov exponent, Ai, Eckmann and Ruelle [7] have 
suggested a certain minimal value for the length of the experimental data points sequence 
in order for the results obtained by applying nonlinear dynamical analysis be reliable. It is 
recognized by several authors that the application of nonlinear dynamical methods implies 
a high computational cost. In this sense, therefore, any method that allows to obtain any 
information about the behavior of a dynamical system and that can operate over a relatively 
short discrete-time series, acquires an important practical value. 
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For complementing those methods of nonhnear dynamical analysis we shall now try to 
demonstrate that introducing the concept of magnitude-squared coherence is possible to 
quantifying the order-chaos transition in a dynamical system described by logistic map- 
ping dynamics. This method is based on the discrete Fourier transform representation of 
a finite-duration discrete-time series and has the advantage that can be applied to a rela- 
tively shorter discrete-time series than that required by nonlinear dynamical methods. The 
magnitude-squared coherence has been used by several authors for measuring constancy of 
phase between two or more signals at one frequency [23,18,24,8,20,16]. Besides, according 
to literature the transition from a regular to a chaotic behavior has been studied in many 
dynamical systems, including the biological type. One of these is a fluid flow, where by 
calculating nonlinear dynamical parameters the order-chaos transition can be determined. 

Van den Bleek and Schouten [4] determined the order-chaos transition in a fluidized bed 
by calculating both correlation dimension, D2, and Kolmogorov entropy, as a function 
of the Reynolds number and the superflcial gas velocity, indicating two different regimes, 
order and chaos, and the intermediate region. In a previous work [9] the DFT representa- 
tion for qualitatively describing the order-chaos transition in a dynamical system is used. 
Now we want to make the description of order-chaos transition in a dynamical system gov- 
erned by the logistic equation dynamics by applying the Welch average periodogram method 
for calculating the power spectral density and from this, to introduce an inverse average 
magnitude-squared coherence index for indicating quantitatively order-chaos transition in 
the dynamical system indicated. 

The application of magnitude-squared coherence (MSG) method is reported by several 
authors. NV Thakor et al [23] presented the power spectral density analysis of electrocar- 
diogram (EGG) waveforms as well as isolated QRS complexes and episodes of noise and 
artifacts. The power spectral analysis showed that the QRS complex could be separated 

4 



from other interfering signals. A bandpass filter that maximizes the signal-to-noise ratio 
would be of use in many monitoring instruments. They calculated the coherence function 
and, from that, the signal-to-noise ratio. Upon carrying out this analysis on experimentally 
obtained ECG data, they observed that a bandpass filter with a center frequency of 17 Hz 
and a Q factor of 5 yields the best signal-to-noise ratio. They estimated the coherence 
function as given by 

(j2 ^ \Gxy{f)-G.xy{f)*\'^ 

Gxx{f)-Gyy{f) 

where G^xif) is the estimates of the auto-power spectra of the QRS complex, Gyy{f) is 
the complete ECG cycle, and Gxy{f) is the cross power spectrum estimate. For calculating 
the power spectra they applied the Fourier transform of the windowing signal. To improve 
the resolution of the FFT algorithm and the power spectra, they filled the original data 
sets with zeros and formed new 1024 point data sets. S Narayanaswamy et al [18] devel- 
oped three signal processing tools to identify three postulated mechanisms of isolated or 
single premature ventricular contraction (PVC) generation. Two hours of continuous ECG 
recording were digitally obtained from several patients with frequent PVC. They applied the 
magnitude-squared coherence (MSG) spectrum between the sinus intervals x[n\ and PVC 
intervals as given by 

MSC{f) = \^-^f^-^y^fy\ (2) 

where Sx{f) and Sy{f) are the Fourier spectra of x[n\ and y[n] respectively. Sy{f)* is 
the complex conjugate of Sy{f), and Sxx{f) and Syy{f) are the auto-power spectra of x[n\ 
and y[n] respectively. The power spectral density of the PVC interval series is defined as 

P{k)^\X{k)\' (3) 

where X{k) is the Fourier transform of x[n]. In this work, details of the used method 
are not showed, but they used both the power spectral density of the PVC signal and the 
MSG estimates of the PVC and sinus intervals. 



5 



EG Lovett y KM Ropela [17] showed that magnitude-squared coherence (MSG) of in- 
tracardiac electrogram may be used to quantify what may be termed organization. Unfortu- 
nately, MSG computations are not meaningful from the surface due to the high correlation 
among measured body surface potentials. They showed that average magnitude-squared 
auto-bi-coherence (AMSABG) provides an organizational measure of cardiac rhythms from 
a single surface lead. They showed that AMSABG also discriminates three ventricular tach- 
yarrthythmia. In this work, rather than measuring consistency of phase between two signals 
at one frequency (as in MSG), measuring consistency between two frequencies in one sig- 
nal (as in MSABG) were considered. Haris J. Sih et al [20] applied the MSG function for 
analyzing the spatial organization of epicardial mapping in the frequency domain. They de- 
fined the magnitude-squared coherence (MSG) as a frequency domain measure of the phase 
consistency between two signals, x[n] y y[n], as given by 

where Sxy{f) is the cross-power spectrum, and Sxx{f) and Syy{f) are the respective auto- 
power spectra, MSG is dimensionless and can vary from to 1. A MSG value of 1 at some 
frequency would indicate a linear relationship or perfect coherence between the two signals 
at that frequency, while a value of would indicate no relationship or perfect incoherence 
between the two signals at that frequency. MSG is sensitive to noise and interference, whereas 
uncorrelated noise in either of the two signals decreases MSG, and correlated interference 
increases MSG. While the relationship between MSG and correlated interference is difficult 
to characterize, the influence of noise on MSG can be derived for two real signals, u{t) 
and v{t), equal to ideal signals, x{t) and y{t), plus random, independent, zero-mean noise 
processes, ni{t) and n2{t). It is easy to show that 

^ SNR:,(f) SNRyif) SNRM)SNRy{f) 

where MSCuv is the magnitude-squared coherence between u(t) and v{t), MSCxy is the 
magnitude-squared coherence between x{t) and y{t), SNRx{f) is the frequency-dependent 
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signal-to-noise ratio between x{t) and ni{t), and SNRy{f) is the frequency-dependent signal- 
to- noise ratio between y{t) and 77-2 (i). As expected, as the signal-to- noise ratios increase, 
MSCuv approaches to MSCxy MSC can be estimated by replacing Sxy{f), Sxy{f), and 
Sxyif) with respective estimates by using the Carter method. This method uses average 
Fast Fourier Transform of weighted overlapped segments of equal length from the two signals 
to form the spectral estimates. That is, MSC is estimated as given by 



where Xi{f) and Yi{f) are the FFT of the i — th weighted segments from x{t) and y{t), 
respectively, N is the total number of segments, and (*) denoted the complex conjugate. 
From this method of estimation, one can interpret MSC as a measure of the phase con- 
sistency between the two signals. At a single frequency /i, the cross-power term in the 
numerator for one segment can be interpreted as a vector with amplitude 
and phase /.Xi(fi) — lYi(fi). Therefore, the sum of the cross-power terms in the numerator 
of the MSC estimate is equivalent to the vector sum of the cross-power vectors for each 
segment. If ZXj(/i) — lYi{fi) — lXj{fi) — Zl^(/i) for all i and j, that is, if a constant phase 
relationship from one segment to next exists between x{t) and y{t) at /i, the N vector terms 
are all aligned and the vector sum is at a maximum, as is the MSC estimate. If the terms are 
not constant from segment to segment, the N vector terms do not align, the vector sum in 
the numerator is not at a maximum, and the MSC estimate must be less than unity at that 
frequency. As with all methods of spectrum estimation, the MSC estimate has statistical 
bias and variance, which have been explicitly calculated for the case of Gaussian stationary 
signals. Another form of bias, known as a bias due to misalignment, has also been shown to 
exist for this estimator. As an example, for two signals, x{t) and y(t) = b* x{t — D) + n{t), 
where n{t) is uncorrelated noise, it has been demonstrated that the estimate of the true 
MSC spectrum will be degraded by a factor of (1 — \D\/Ty where T is the time duration 
of one windowed segment of data. For example, if T were approximately 0.25s, a relative 
delay of 0.1s between the two signals could decrease the MSC estimate to one-third the 
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true value. It is noted that these results are derived for two linearly related signals, and 
the effects of this bias due to misalignment on two nonlinearly related signals are yet to be 
determined. With various types of mapping data, as those of cardiology, of neurology, of 
geophysics, etc, multiple signals are recorded simultaneously, and often, it is necessary to 
preserve the physical location of these recordings relative to the underlying activity. Under 
this conditions, n signals, xi{t),X2{t), -jXnit), are arranged in a 1—, 2—, or D— array, or 
map. By choosing some meaningful reference signal, Xr{t), which can be one of the n signals 
in the map, a corresponding map of MSG spectra can be formed. That is, what we have is 
an — MSC spectra set: MSCxj^xrif)^ ^SCx2Xr{f)-i' " ") MSCx„xr{f)j arranged in the same 
1— , 2—, or D— map. In the case where the reference is one of the n signals in the map, the 
r — th MSC spectrum, MSCx^xXf)-, unity for all frequencies. It is noted that the map 
of MSC spectra may indicate how the n signals relate to the reference, but they do not 
imply any relationship between one signal of the n signals with any other of the n signals. 
Different maps of MSC spectra can be generated using different references. A logical and 
concise method for displaying these MSC spectra is problematic as pointed out by many 
authors. To simplify the data, one can reduce each MSC spectrum to a single number, for 
example, by averaging MSC over frequency. 

Finally, we refer to the work done by D Liberati et al [16] , where the joint use of total and 
partial coherence between pairs of EEGs simultaneously recorded in a standard set is shown 
to enhance what is caused direct correlation between cortical subsystems and what is instead 
related to the spread of the electromagnetic field. A multi- variable AR approach is employed 
in the computation, giving results even for a very short time window, thus allowing coherence 
to be investigated at the main cortical latencies of evoked potentials. In particular, when a 
combined visual and somatosensory stimulation is applied, cortical interactions are captured 
in the frequency domain. In this work, the total squared coherence A;, as a function of the 
frequency / is defined, for every couple j, h of the channels, as the squared cross-spectrum 
Pjhi normalized to the product of the norm of each auto-spectrum Pjj and Phh as given by 
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where kjh is a real function, ranging from to 1, indicating the amount of correlation 
between the signals j and /i as a function of the frequency /. In this work we pretend to 
demonstrate that it is possible to apply the average Welch periodogram for estimating the 
power spectral density of a chaotic signal and to calculate the magnitude-squared coherence 
for determining quantitatively the order-chaos transition in a dynamical system described 
by the logistic mapping dynamics. The fundamental advantage of this method upon the 
others is that it can operate on a relatively shorter discrete-time series. By estimating the 
largest Lyapounov exponent was able to evaluate those results obtained here applying an 
inverse average magnitude-squared coherence index, defined by averaging MSG estimate 
over the frequency. 

II. METHODS AND MATERIALS 

A. Obtaining the discrete-time series set 

The logistic equation is an one-dimensional mapping of the real axis in the interval (0,1) 
to itself, that is a prototype of a nonlinear dynamical system widely used [6] and can be 
formulated as 

-|- 1] = ra;[n](l — (8) 
and it can also be formulated as 

n 

x[n + 1] = r"(x[l] - J2 r~''^^x[k + 1]); (9) 

fc=0 

and with the initial condition x[l] specified, one can obtained a time sequence of a given 
length by a recursive action. The character of the behavior of a data sequence obtained 
with the iterative process can be modified by conveniently selecting a value of the control 
parameter r. In this work a value of the initial condition x[l] — 0.65 was taken, and the 

9 



interval r e [2.8, 4.0] was selected. A set of evenly spaced values of r, with a step S — 0.01, 
was taken to produce a family of 120 discrete-time series, each containing N — 2^° samples. 
The building up of the histogram for each time series allows to have immediately a simple 
statistical characterization of this time series. In order to make easier further computations 
data points were organized into a rectangular matrix containing 120 columns, each of which 
is a time series with the length N = 1024 data points. 

B. Determining the largest Lyapounov exponent Ai 

Lyapounov exponents spectrum deals with average exponential divergence or convergence 
of two neighbor orbits in the phase space. It is common to order the spectrum from largest to 
shortest exponents: Ai, A2, A3, .... A system containing one or more positive exponents is to 
be defined as chaotic. For calculating the Lyapounov exponents spectrum some algorithms 
have been developed, according to the particular situation, being one of the most used 
algorithms that reported by Wolf et al [25]. Prom the exponents spectrum, the largest 
exponent, Ai, decides the behavior of the dynamical system. For calculating the largest 
exponent the algorithm proposed by Rosenstein et al [19] can be referred too. Details of 
these algorithm can be found in referred papers. We determined the largest Lyapounov 
exponent, Al, for each time series in the set, using a professional software [21]. The selection 
of the time delay, r, and the embedding dimension, D^, required for reconstructing the phase 
space is part of the problem and it is necessary to consider some selecting criteria. This 
selection can be made by different methods [4,22,13,10,7,15,11]. For selecting r we may 
take into account: (a) optimal filling of the phase space [5], (b) by determining the position 
of the first local minimum of the data autocorrelation function [1], (c) by taking the first 
minimum in the mean mutual information plot [12], (d) based on the optimal tradeoff 
between redundance and the irrelevance [19], and by analyzing the differential equations of 
the system if they are given [26]. For selecting De it may be indicated: (a) by applying 
the Grassberger-Procaccia method, which allows to obtain simultaneously both correlation 
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dimension and embedding dimension [13], (b) by applying the false nearest neighbors method 
[19,3], and (c) by applying geometric considerations [15]. In this work, a time series obtained 
from the logistic equation was analyzed by using the method (a) for r determination and 
by taking a value D^. — 3, which satisfies the Grassberger-Procaccia algoritm condition 
Df. > 2D2, and Takens theory [13,22], being D2 — 0.5 a value reported by Hoyer et al [14] 
for this time series type. Literature refers [13,22] that for selecting the embedding dimension 
as 

De>D2 + 1, (10) 

it is sufficient to have a good representation for the attractor in the phase space, re- 
constructed in delayed coordinates from a scalar time series observed. Any reconstruction 
made with a value of the embedding dimension less than the minimal value will produce a 
projection, in that dimension, of the original attractor, and therefore it will be more diffi- 
cult to understand that projection or it will produce false results when an estimation of a 
nonlinear dynamical parameter is to be made. It is important to take a value for embedding 
dimension close to 2D2 -|- 1 in order to avoid the computational cost be excessively high. 

C. Estimating the lAMSCI using the Welch average periodogram 

It is highly meaningful when dealing with a chaotic system keeping in mind that lo- 
cal sensitivity to a small error is the hallmark of such a system. Such dynamical system 
as the one described by the logistic equation produces a discrete-time non-chaotic series 
when the control parameter r value satisfies the condition r < 3.57. For such a condition a 
discrete-time series generated by the recursive process does not almost change if the initial 
condition is modified and there exists phase consistency between the two discrete-time series 
at a given frequency. However, when the control parameter value is greater than 3.57, a 
discrete-time series generated by the logistic equation does meaningfully change if the initial 
condition is slightly modified and phase consistency between the two discrete-time series at 
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a given frequency, in general, must be lost. This effect may be quantified by applying an 
inverse average magnitude-squared coherence index. In order to apply the Welch average 
periodogram method for determining the power density spectrum Sx[k] of a discrete-time 
series one follows the following steps: 



(a) Decompose the sequence of N (=1024) data points in M (=9) segments, each one 
with the same length L (=256), and which may be overlapped P samples. The usual value 
P — 0.625 * L was chosen. 



(b) Calculate the DFT representation by FFT algorithm for each segment as given by 
-'^mM = fftixjnln]), where m = 1, 2, . . . , M, and n = 0, 1, 2, . . . , L — 1, and fft represents a 
function in MatLab for calculating the DFT of a given signal. For each segment a L-FFT 
was calculated, being L the number of points in the sequence Xm[^]- 

(c) Calculate the periodogram for each segment m as given by 

Pm[k] = ^\Xm[k]\\ m = l,2,...,M (11) 

(d) Finally, calculate the average periodogram as given by 

1 M 

-mT. Prn[k\ (12) 
-''^ m=l 

where k = 0, 1, 2, . . . , L — 1. For a short discrete-time series a rectangular window is 
recommended. For determining the magnitude-squared coherence sequence the procedure 
described above is applied to a second signal y[n] obtaining the average periodogram Sy[k], 
and the average cross-periodogram S^ylk] is also calculated. The magnitude-squared coher- 
ence sequence was calculated as given by 

The MSCxy[k] sequence describes the phase consistency between the two signals at a 
given frequency k. The mean value of the sequence, {{MSCxy[k])), may be used as an index 
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of global coherence between the two signals. When deahng with a chaotic discrete-time 
series this index must tend to zero, and when dealing with a deterministic discrete-time 
series this index must tend toward one. However, it was found that a better index may be 
defined as given by 

It makes easier to compare the results obtained here with those obtained by applying 
the largest Lyapounov exponent estimate. 



III. RESULTS 

Figures, from 1 to 4, show a sample of four representative time series and their cor- 
responding magnitude-squared coherence spectra. In order to give a better view of each 
time series it was only considered 128 data points of each discrete-time series in its plotting, 
but it was completely considered when its corresponding magnitude-squared coherence was 
calculated. In either case, the particular values of r parameter were indicated. Figure 5 de- 
picts results of plotting the inverse average magnitude-squared coherence {1/ {MSCxy[k])), 
calculated for each time series in the set, as a function of the control parameter r. And 
figure 6 shows the results of plotting the inverse average magnitude-squared coherence in- 
dex, in decibels (/ = 10 * logio{l/ {MSCxy[k])), calculated for each time series in the set, 
as a function of the control parameter r too. Figure 7 shows results of plotting the largest 
Lyapounov exponent estimate, Ai, calculated for each time series in the set, as a function of 
the control parameter r. This reference serves as a control for the discussion of the results 
obtained by applying the inverse average magnitude-squared coherence index method. 



IV. DISCUSSION OF THE RESULTS 



The application of the DFT representation of a chaotic-type signal gives satisfactory re- 
sults for determining the order-chaos transition in a system described by the logistic equation 
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dynamics [9] . A critical value for control parameter r is reported in literature beyond which 
a sequence produced by the logistic equation exhibits a chaotic behavior. This threshold 
value approaches r = 3.57. On the other hand, figures from 1 to 4 show the time domain 
representation of some of the time series and their corresponding magnitude-squared co- 
herence spectrum. Comparing figures 1 (r = 2.91) with 2 (r — 3.64), 3 {r — 3.75), and 
4 {r — 3.84), it can be deduced the bifurcation effect and the transition to chaos. Note 
that in figures 2, 3 and 4 the control parameter satisfies the condition r > 3.57, and it 
corresponds to a value for a positive largest Lyapounov exponent, as it can be estimated 
from figure 7. It may be observed in those cases that the magnitude-squared coherence 
spectrum exhibits values far away from one for several values of the frequency in the MSG 
spectrum. It makes the average magnitude-squared coherence to be less than one, corre- 
sponding to a situation for which the discrete-time series becomes incoherent and indicating 
in this case a chaotic behavior of the dynamical system. This situation can be observed 
in figures 5 and 6 where either the quantity {1/ {MSCxy[k])) as the inverse average MSG 
index, (/ = 10 * logio{l/ {MSCxy[k])), are greater than 1 and 0, respectively, for r > 3.57. 
Applying the inverse average magnitude-squared coherence index (lAMSCI) to a discrete- 
time series in a family derived from an observable in a dynamical system is recommended 
for quantitative detecting the order-chaos transition, and this can be added to the metric 
tools of the nonlinear dynamical analysis for complementing. Any method that allows to 
evaluate the behavior of a dynamical system, and that can be applied to a relatively shorter 
discrete-time series, must be taken into account as an important practical method. 



V. CONCLUSIONS 

The apphcation of the inverse average magnitude-squared coherence index (lAMSCI), 
/, to a discrete-time series obtained from the logistic equation can be done on a chaotic 
discrete-time series. This discrete-time series has the property that exhibits a high level of 
loss of coherence, and it can be detected by applying the lAMSCI method. This method 



14 



does not exert a strong requirement on the length of the experimental data sequence. If 
one deals with a nonlinear dynamical system which can modify its behavior between order 
and chaos, a time series produced by that system reflects this change, and the lAMSCI 
method applying to the discrete-time series allow to determine this change. Any method 
that permits some quantitative evaluation about the behavior of a dynamical system and 
that can work over a relatively short time series acquires a practical importance. Besides, 
this method does not require the phase space reconstruction of the attractor. 
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FIG. 1. Discrete-time series and its magnitude-squared coherence, for r=2.91 
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FIG. 2. Discrete-time series and its magnitude-squared coherence, for r=3.64 
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FIG. 3. Discrete-time series and its magnitude-squared coherence, for r=3.75 
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FIG. 4. Discrete-time series and its magnitude-squared coherence, for r=3.84 
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FIG. 5. Inverse average magnitude-squared coherence index 
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FIG. 6. Inverse average magnitude-squared coherence index, in dB 



22 



X Experimental value 
Behaviour 



Control parameter, r 



FIG. 7. Largest Lyapounov versus control parameter r 
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